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Introduction 

Many applications require an algorithm that averages quaternions in an optimal manner. 
For example, when combining the quaternion outputs of multiple star trackers having this 
output capability, it is desirable to properly average the quaternions without recomputing 
the attitude from the the raw star tracker data. Other applications requiring some sort 
of optimal quaternion averaging include particle filtering [1] and multiple-model adaptive 
estimation [2], where weighted quaternions are used to determine the quaternion estimate. 

For spacecraft attitude estimation applications, [1] derives an optimal averaging scheme 
to compute the average of a set of weighted attitude matrices using the singular value de- 
composition method [3]. Focusing on a 4-dimensional quaternion Gaussian distribution on 
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the unit hypersphere, [4] provides an approach to computing the average quaternion by min- 
imizing a quaternion cost function that is equivalent to the attitude matrix cost function 
in [1]. Motivated by [1] and extending its results, this Note derives an algorithm that deter- 
mines an optimal average quaternion from a set of scalar- or matrix-weighted quaternions. 
Furthermore, a sufficient condition for the uniqueness of the average quaternion, and the 
equivalence of the minimization problem, stated herein, to maximum likelihood estimation, 
are shown. 

For the scalar weighted case the goal is to find the average of a set of n attitude estimates, 
q, , in quaternion form with associated weights W{. The simple procedure 

( n \ ~~ 1 n 

Y^ w i) Yl Wiqi co 

i= 1 / i = 1 

has two flaws. The first and most obvious flaw, that q is not a unit quaternion, is easily 
fixed by the ad hoc procedure of dividing q by its norm. The second flaw is subtler. It is 
well known that q and — q represent the same rotation, so that the quaternions provide a 
2:1 mapping of the rotation group [5]. Thus changing the sign of any q; should not change 
the average, but it is clear that Eq. (1) does not have this property. 

The observation that we really want to average attitudes rather than quaternions, first 
presented in [1], provides a way to avoid both of these flaws. Following this observation, 
the average quaternion should minimize a weighted sum of the squared Frobenius norms of 
attitude matrix differences: 


( 2 ) 


q = argmin VV \\A(qj - A(q<)||J 

q €§ 3 Z — s 

i—1 

where S 3 denotes the unit 3-sphere. 

The Average Quaternion 


Using the definition of the Frobenius norm, the orthogonality of kl(q) and A(q,), and 
some properties of the matrix trace (denoted by Tr) gives 


IKq) - A(cn)f F = Tr |[A(q) - M( qi )] T [A{ q) - ^(q*)]} 
= 6 -2 TV [kL(q)A r ( qi )] 


( 3 ) 
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This allows us to express Eq. (2) as 


q = argmaxTr [A(q) B T ] (4) 

where 

n 

B = ^2wiA(qi) (5) 

i= 1 

Equation (4) is in a form found in solving Wahba’s Problem [6] , so many of the techniques 
used for solving that problem [7] can be applied to finding the average quaternion. If compu- 
tational efficiency is important, the well-known QUEST algorithm [8] can be recommended, 
as will be discussed later. The matrix B is known as the attitude profile matrix [9] since it 
contains all the information on the attitude. 

A detailed review of quaternions can be found in [5], but we only need a few results for 
this paper. We denote the vector and scalar parts of a quaternion by q = [g T <74] T , which 
are assumed to obey the normalization condition ||^|| 2 + g| = q T q = 1. The attitude matrix 
is related to the quaternion by 


X(q) = (q\ - Hell 2 ) / 3 X3 + 2e e T - 2<j4|ex] 


( 6 ) 


where / 3x3 is a 3 x 3 identity matrix and [px] is the cross-product matrix defined by 


[f?x] 


0 -<73 <72 

<73 0 — <7i 

-<72 <7i 0 


(7) 


Equation (6) can be used to verify the identity 


Tr[A(q) B T ] = q T iPq 


( 8 ) 


where K is the symmetric traceless 4x4 matrix 


K 


B + B T - Tr{B)I 3x 3 z 


Tt(B) 


(9) 


with z being defined by 


zx]=B T -B 


( 10 ) 


This is the basis of Davenport’s q- method [7]. The case at hand admits considerable simpli- 
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fication, however. Substituting Eq. (6) for A(qj) into Eq. (5) and then into Eq. (9) gives 


K — 4 M Wtoi^AyA 


( 11 ) 


where w tot = Xu=i -^ 4 x 4 is a 4 x 4 identity matrix and M is the 4x4 matrix 

M -5Z w *q*qr (12) 

Thus the average quaternion can be found by the following maximization procedure: 

q = argmaxq T Mq (13) 

qeS 3 

The solution of this maximization problem is well known [10]. The average quaternion is 
the eigenvector of M corresponding to the maximum eigenvalue. This avoids both of the 
flaws of Eq. (1). The eigenvector is chosen to have unit norm to avoid the first flaw. The 
second flaw is obviously avoided because changing the sign of any q; does not change the 
value of M. The averaging procedure only determines q up to a sign, which is consistent 
with the 2:1 nature of the quaternion representation. If the QUEST algorithm is used to 
find the eigenvector associated with the maximum eigenvalue, then the K matrix in Eq. (11) 
must be used instead of the M matrix in Eq. (12), because the QUEST algorithm requires 
a traceless matrix. 

A closer look at the attitude error matrix defined by A((5qQ = A(q)A T (qj) gives a nice 
interpretation of the optimization problem. The error quaternion is the product of q and 
the inverse of q*, which can be written as [5] 

<5q; = q <g> q" 1 = [H(qQ q;f q 

where 

S(q) = 

Note for future reference that [H(q,) q*] is an orthogonal matrix representing a norm- 
preserving rotation in quaternion space. The vector and scalar parts of the error quaternion 
are given by 


94^3x3 + [£ x ] 


(15) 


(14) 


S@i = sm(S(f>i/2) = E T (q Q q 

5q M = cos{5(j)i/2) = qfq 


(16a) 

(16b) 
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with e,- being the unit Euler axis and <50; the rotation angle of the error. Substituting 
Eqs. (16) into Eqs. (2) and (3) and using the quaternion normalization condition shows that 
the average quaternion is given by 


q = arg min ) W;||<5£;|| 2 =argmin> sin 2 (<50;/2) 

< J q €§ 3 1 4 


qeS 3 . 
1=1 


The fact that the average quaternion minimizes the weighted sum of the squared lengths of 
the vector parts of the error quaternions, or the weighted sum of the squares of the sines 
of the half-error-angles, may be more intuitively pleasing than the argument based on the 
Frobenius norm. Equation (17) will be used in the sequel for the generalization of the results 
to non-scalar weights. 

Averaging Two Quaternions 

The two quaternion case can be solved in closed form. The optimal quaternion average is 
given by 


q = ± 


[(u>i - W2 + *)qi + 2w 2 (q r fq 2 )q 2 ] 
\\(wx-w 2 + z) qi + 2 tu 2 (qfq 2 )q 2 | 
[2wi(qfq 2 )qi + (w 2 - Wi + z)q^ 
||2wi(q'fq 2 )qi + (w 2 -Wi+ z) q 2 | 


where z = \J (w[ — w 2 ) 2 + 4witt) 2 (qfq 2 ) 2 . These two forms are equivalent if qfq 2 ^ 0. If 
qfq 2 = 0, then z =\w\ — vj 2 j ; the first form gives the correct limit q = qi if w\ > w 2 , while 
the second form gives the correct limit q = q 2 if W\ < w 2 . If qf q 2 = 0 and 'u>\ = w 2 , neither 
form has a well-defined limit. This is true because the maximum eigenvalue of M, which is 
equal to wi — w 2 , is not unique in this case. This is the only two-observation case for which 
the average quaternion is not uniquely defined. 


Uniqueness of the Average Quaternion 

Because the average quaternion q is the eigenvector associated with the maximum eigen- 
value of M, q is unique if and only if the two largest eigenvalues of M are not equal. A 
sufficient condition for the uniqueness of the average quaternion is shown here. It is assumed 
there is a reference frame in which every quaternion estimate q; differs from the identity 
quaternion q re f =[0 0 0 1] T by a rotation of less than tt/2. This section proves that the 
average quaternion q minimizing Eq. (2) is unique with this assumption. 

The angle of rotation between q; and q re f is given by 2arccos |g 4 J. When the angle is less 
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than 7r/2, q\. > 1/2, hence 


( 19 ) 


2 ^ 2 ! 2 . 2 
?4, > 9i, + 92, + 9 3< 

Now consider an attitude quaternion that is orthogonal to the identity quaternion q re f, given 
by q- 1 = [q ± q% q$ 0] T . Define the gain function g( q) = q T Mq. The gain functions of q- 1 
and q re f are 


n 

w i(Qi Qii + 4 (hi + 4 hi) 2 (20) 

i= 1 


and 


5(qref) = ^ Wi 4 ( 21 ) 

i — 1 

We have 

9 (qref) > g{ q 4 (22) 

because 

4 > (4 qu + 4 q*. + 4 q^f ( 23 ) 

To prove inequality (23), notice that, by the Cauchy inequality, 


(4 qu + 4 9 . 2 , + 4 q^) 2 < [(4) 2 + (4) 2 + (4) 2 ] (4 + 4 + 4) 


2 . 2,2 

q u + q 2i + q^i 


(24) 


where the fact that q x has unit norm has been used. Inequality (23) then follows upon 
combining Eqs. (19) and (24). 

Now, if the two largest eigenvalues of M are equal, the eigenvectors associated with the 
maximum eigenvalue span a 2D subspace. The intersection of this subspace and the orthog- 
onal complement of the identity quaternion (the subspace spanned by all q x quaternions) 
cannot be empty. A quaternion q in that intersection must satisfy y(q) = </(q) > f/(q re f) 
and y(q) < g(q re f) simultaneously. By contradiction, the two largest eigenvalues cannot be 
equal, hence the average quaternion is unique. 


Matrix Weighted Case 


This section expands upon the scalar weighted case to include general matrix weights. 
For this case, a matrix weighted version of the minimization problem in Eq. (17) is assumed: 


arg min yw 
i = 1 


i R i 


1 6qJ = arg min q T 5(q,.) R t X S T (q *) q 

qG§ 3 ' 
i—l 


(25) 
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where R { 1 is the z th symmetric weighting matrix. In this case the average quaternion is the 
eigenvector corresponding to the maximum eigenvalue of the matrix 

m = - ( 26 ) 

i= 1 


If Ri 1 = Wil 3 X3 , the identity S(qj) ^ T (qi) = UvA — q 2 qf can be used to show that Ri = 
M — w to t hx4i where M is given by Eq. (12). The traceless K matrix for the matrix- weighted 
case is 

RT 1 1 hxi = 

) i=i 


£ K < ( 2? ) 


K 




, i = 1 


with 

Ki ± -4 E( qi ) i?r 1 E T (q i ) + Tr (R; 1 ) I 4x 4 = [H(q») q i] Ki [H(q 0 q^ (28) 


where 


2 Ti~ Tr(jFj)/ 3X 3 


Ki = 2 


0 3 xl 


OL, 


(29) 


(30) 

and 0 3xl denotes a 3 x 1 vector of zeros. 

The matrix Ki has the same structure as the one in Eq. (9) with the corresponding 
attitude profile matrix given by Bi = 2T{. Because [S(qj) q,] is an orthogonal matrix, 
Eq. (28) shows that K % and K{ are related by a rotation. Their corresponding attitude 
profile matrices are related by the same rotation. The attitude profile matrix corresponding 
to Ki is given by Bi = Bi A(qi), as is shown in the Appendix, and the attitude profile matrix 
corresponding to K is given by 


b = E b ‘ = 


2E^(q.) 

4=1 


(31) 


If Ri 1 = w t I 3x3 , Eqs. (27) and (31) reduce to the corresponding quantities for scalar weights. 


Relation to Maximum Likelihood Estimation 

The relationship of the minimization problem in Eq. (25) to a maximum likelihood estimation 
problem is now shown. Reference [11] establishes a maximum likelihood problem for attitude 
matrices, which is related to the averaging problem in this paper. The maximum likelihood 
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estimate of the attitude matrix, denoted A u l, is given as 


i A .1 

A M l = arg nun - 
Aeso 3 2 


J^TrKAi-AfFiiAi-A)] 


i — 1 


(32) 


where Ai is the z th given attitude matrix, SO 3 denotes the (special orthogonal) group of 
rotational matrices and the matrix T % defined by Eq. (30) is the Fisher information matrix 
of the small attitude matrix errors, with R t being the covariance of the small attitude vector 
errors. Using properties of the matrix trace, we can write 

Ji(A) 4 I Tr [(Ai - A) T F t (Ai - A)] = Tr^) _ I Tr (A Bj) (33) 

where Bi is the attitude profile matrix defined by Eq. (31). Using A = A(q) and the definition 
of J~i gives Ji as a function of the quaternion: 

J,{q) = Tr(^) - i Tr p(q) Bj] = i [Trfi?" 1 ) - q T K, q] = 2 SefRr'SeJ (34) 


where Eqs. (8) and (28) have been used. Using the invariance property of the maximum 
likelihood estimate [12], 

^4-ml = 2f(qMiJ (35) 


where qML is the maximum likelihood estimate of the quaternion. Hence, using Eq. (34) in 
Eq. (32) gives 


qML = axgmm^ 5 gjRi Agf 




(36) 


which is identical with Eq. (25). Thus, we conclude that the average quaternion defined by 
Eq. (25) is a maximum likelihood estimate. 

The error-covariance associated with the small-angle attitude errors of the average quater- 
nion is given by 


R 


t:T 


(q) 


^Eiq^Ri 1 E T (q i ) 


i=l 


-1 


=(q) 


(37) 


For small errors, this matrix is well approximated by 




(38) 


, i=l 


Equation (38) can be used to develop 3-sigma bounds for the attitude errors between the 
average and true quaternion. 
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Conclusions 


An algorithm is presented for determining the average norm-preserving quaternion from 
a set of weighted quaternions. The solution involves performing an eigenvalue/eigenvector 
decomposition of a matrix composed of the given quaternions and weights. For both the 
scalar- and matrix-weighted cases, the optimal average quaternion can be determined by the 
computationally efficient QUEST algorithm. A uniqueness sufficient condition is presented 
for the scalar-weighted case. In the matrix-weighted case, when the matrix weight is given 
by the inverse of the covariance of the small attitude-vector errors, the average quaternion 
is shown to be a maximum likelihood estimate. Thus, in this case, the averaging procedure 
introduced here enjoys the well known desirable properties of maximum likelihood estimators. 

Appendix 

This appendix proves that Bi = B l A(qj). Here we let unsubscribed q denote a normal- 
ized, but otherwise completely arbitrary, quaternion. Then from Eqs. (8) and (28), 


Tr[A(q )Bf] = q T K t q = q T [S(q;) q 4 ] K t [S(qQ q*] 7 ’ q (A-l) 

Using Eq. (14) and the fact that the attitude profile matrix of iQ is Bi gives 

Tr[A(q)Bf] - dqf Ki dq, = Tr[A(6qi)Bf} (A-2) 

It follows from Eq. (14) that A(Sqj) = A(q®qj~ 1 ) = A(q)A T (qj), so 

Tr[A(q )BJ] = Tr[A(q)A T (qi) Bj] (A-3) 

Since this must be true for any quaternion q, it follows that 

B i = [A T (q i )Bn T = B l A(q i ) (A-4) 
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